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1, Introduction. Numerous papers developing the theory of 

interpolation and approximation by spline functions have ap- 

peared in recent years. Concurrent with this theoretical de- 

velopment, there has been considerable interest in practical ; 
algorithms for the computation of splines. Unfortunately, how- ; 
ever, agreat deal of the work on algorithms has not been doc- 

umented, or at best is available only in local reports of the 

academic departments, and industrial, government, or private 

laboratories in which it was carried out (cf. the references 

cited at the end of this article). Thus, a great wealth of com- 

putational experience and accompanying insight is not so read- 

ily available to the prospective user of spline methods. 

The purpose of the present survey is to at least make 
the reader aware of the existence of various algorithms for 
computing interpolating and approximating splines, to outline 
their derivation whenever practicable, and to make some at- 
tempt to discuss the numerical experience obtained with them. 
Because of the above-mentioned relative inaccessibility of 
sources, it is essential that a strong disclaimer of any kind 
of completeness be made at the outset. 

The first part of the paper is devoted to soine typical 
interpolation problems and corresponding algorithms for their 
solutions, while the second part is concerned with approxima- 
tion methods. " Whenever possible, examples of the kinds of 
problems the algorithms were successfully tested on are quoted, 
as well as some indications of where failures occurred. The 
details of programming the algorithms are completely ignored, 
and no existence, uniqueness or convergence results are es- 
tablished in this survey. The reader is strongly advised to 
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consult the available references as well as experienced prac- 
titioners if possible before attempting large scale computa ~ 
tional use of what appears here. 


2. Interpolating splines. 


2.1 Interpolation at nodes differing from the knots. A spline 
function of degree n with k prescribed knots RS Ko Swe Ss xX, 


is a function s(t)e¢ C®-! such that s(t) € mp, the class of 
polynomials of degree at most n , in each of the intervals 
(90,2), (%),%5),+++,(X,,%)- Denote by Sn, k( Xp +++ 5X) 
this class of splines. 


Problem 2.1 Given knots xX) SxX2 <.-.< x, points 
ti <tg<... <thykq) satisfying 


; <i eT ee Toenar 
(2.1) as Lee fet By wate ye 


and real numbers ga +» compute a spline function of 


degree n with knots {x,}§ interpolating {y, pret at the 
 ielaid " 


We first remark that since {l,t,...,t™, (t - 
(t-x,),} form a basis for Sy K(X, ++) XK), Coe. 
Ss must satisfy 


n 
is 


n k 
i at n 
(2. 2) 9(t)) 2 a, t+ 2, elt x) Hy, 
j=1,2,...,ntk4+l. 


Since the determinant of ( 2.2) is nonzero (see Lemma 2. 2 of 
[22]), s exists uniquely. 

The simplest possible algorithm for computing s 
would amount to solving the linear system (2.2) by any of the 
usual methods. This turns out to be impossible in practice, 
however, since the matrix of the system is very badly condi- 
tioned unless n and k are both small. To derive a viable 
algorithm for obtaining s numerically it is necessary to find 
a different basis for Sn, K(X p20 Xk) - 

Choose X-n<X-pn4]<... <xQ < x) and xk <xXp4) < 

++ SXkgn¢1- For i =1,2,...,mtk+l define the B-splines 
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sp aeel 


(2.3) M,(t) = Mts, 1p Xyige ry 


where M[-] denotes the (n + 1) th divided difference on the 
points xj -y-}, Xj-n’***x; Of the function M( tiy)= (n+]) 
(y - t)P. By Theorem 5 of Curry and Schoenberg [6], every 
s(t) € $n, k(x), ..+,Xk) has a unique expansion 

ntk+1 


(2,4) BCU c, M(t) for x) <t<x 


=l te) k+l ° 


This basis for 8, KOXp a x) leads to the following con- 
venient interpolation method. 


Algorithm 2.2 (Compute the coefficients tape in 


the expansion (2.4) for the solution of Problem 2.1). 
(a) Set h =(x,-x)Ak-1), x9 =min(ty,x; -h) , 


Xk4] = Max (tnek+], X +h), and x_j =xQ -ih 
Xiekel = Xk + th for Pes, Boas nahin. 


(b) Solve the system 


ntk+1 
(2.5) 2 cM (t)) =¥,, 


{= 1 Bin opnbler 


for the ie ee 


This particular interpolation problem has been all but 
neglected, and this seems to be the first time an algorithm 
concerned with it has appeared in print. The method is easily 
modified to allow the computation of Tchebycheffian splines 
(see §6 of [22]) interpolating as in Problem 2.1. 

Algorithm 2.2 has been programmed in Algol at the 
Mathematics Research Center for the Burroughs 5500 computer, 
and has been tested on a number of examples. Good results 
for fairly wild data (y;'s oscillating up to say 104) have been 
obtained for modest sizes in the degree of the splines, e.g., 
up tolO. For high degrees some difficulty in accurately com- 
puting the divided differences involved in obtaining the M,(t)) 
was encountered. It should be remarked that the matrix in 
(2.5) is of the "banded" type because of the restricted sup- 
port of the Mj's , and this local property of the expansion 
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is also favorable in evaluating the spline. 


2.2 Natural splines. A function s( t) is called a natural 
spline of degree 2m-l with knots x)<x.<... <x, provided 
Se€ S2m-l,k(%)) seeyX%p) (cf § 2.1) and s(t)e™p,-) in each 
of the intervals (-90, x)) é (x), oo), 


Problem 2.3 Given x; <x2<... <x, and real num- 
bers {yj y compute a natural spore function of degree 2m-l 
interpolating fy} at the {x,}} 

The existence of a unique s satisfying Problem 2.3 
(m<k) follows from the fact that it is determined by the sys- 
tem of equations (cf [15]) 


oa i F am-l 
(2.6) atx.) s a, se 4 Gog Se) =Y,, 
j i=0 oh yee Mea j 
fh, Shc0s eh 
¥, r 
ey cx =0 Pre O ewe Ot 


As in §2.1, the simplest algorithm for solving Problem 
2.3 would involve numerical solution of (2.6), but the same 
objection applies; it is an ill-conditioned system. Several 
alternate approaches have been devised. 





Let n= 2m-1 and choose {x,}%, and {x; fas as 
Xy-j =X, ~ th, xXjy, =Xk + th, i =1,2,...,m4l, where h= 
a ke Let {m,( t) ae be the B-splines defined in 
(2.3), and let (2.4) be the expansion of s. Now 
ntk+1 
(zc?) s(x;) * co; M, (x;) 259 j=1,2,...,k 


provide k equations for the determination of the ey gander 
Consider dividing the intervals [x g,x,] and [Xk, X, ] into 
n equal parts. Since s( m) ( t) =O in these intervals we have 


»m-l 


(2.8) f a™s(x, +) hg) =0 {S05 0 


AM s(x, +i hy) =0 j =0,1,...,m-l 
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where hg =(x,~X9)/n, hy =(x,4) —xK)/n, 

and the differences are taken with respect to j. The equa- 
tions (2.7) combined with (2.8) provide n+tk+l equations 
for the ntk+l c's. 


Algorithm 2.4 (Determine coefficients ii in 
the expansion (2.4) of s solving Problem 2, 3). 
fetes x - ih and Xiay =X) + ih 
for (fa, 2). RN. Th Sie ig) CK = Ue 


b) Compute M,(x;) 5 EE 25 5 eo ghey, Sle Zine eat es 


a) Set x 


c) Compute AM s(Xq + j ho) 7 20, 1.2. pak; 


A™ s(x, + j hg) hg =h/n. 

d) Solve system (2.7), (2.8) for inne i 

Algorithm 2.4 has also been tested at the Mathematics 
Research Center on a number of examples, again with fairly 
wild data. As with Algorithm 2.2 good results were obtdined 
modulo some difficulty in computing the divided differences 
for higher degree splines. The matrix here is also fairly 
sparse, and the local property of the basis likewise obtains. 

Still another method for computing the solution to Prob- 
lem 2.3 relies on obtaining an expression for sl m)(t), We 
need to introduce the notation 





PNP 5+ 19 My _ a ; 
ae ce a iz ( 25-2). (2j-Z,_y) (4-25) (2-2) 


for the m-th divided difference of the data 9), Pay 06s Pm at 
the points 2),22,...,2m. It can be shown that 


=] % - 
k=m dM; (t) m-l t 


(2.9 s(t) = es evaee e eG t 
) (%) ia. (amt * Yo 4% 
where 
2m-l 2m-] 
(x,t) rece. (i507 


(2.10 ) M, (1) = mD 
Ey 
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Algorithm 2.5 (Determine ior? ! igie™ in the 
expression (2.9 ) for the solution of Problem 2.3.) 


Yyrre eo yan 

497 9 am 
Evaluate My (x) 5 Gel Aiding ky Jely aia s earn: 
Evaluate Ay by (2.13) 4,.f2),.2,. 5.0, Kem . 

b) Solve (2.12) for ae 

c) Solve (2.11) for {q,)a-! ‘ 


a) Evaluate D 





| fi, 2,..05.,.Rom: 


It may be noted that step c) amounts to determination 
of the unique interpolating polynomial of degree m~1 for m 
given data points, 

This algorithm has been tested on the same data as 
Algorithm 2.4, with comparable results, although because of 
the nonlocal character of the underlying basis, some loss in 
accuracy near the right-hand end of the interval was noticed 
for large numbers of knots. Moreover, accuracy was found to 
suffer when extremely nonuniform knot spacings were pre- 
scribed, e.g., with the ratio of largest to smallest knot dis- 
tance on the order of 100 or 1000. 

Numerous authors have obtained Algorithm 2. By ior 
slight modifications of it, see, e.g., [1, 2, 3, 4, 5, 6, ll, 
14, 16] etc. The above form was taken from [16], where it 
is shown to also be applicable for Tchebycheffian splines 
(see §6 of [22]). Considerably more general splines 
have been studied, and analogous algorithms for their 
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computation are available, see [17] and references 
therein, 

The idea of a smoothing spline (see [18 , 19 ]) is 
very closely related to interpolating splines, and there is an 
algorithm in the literature [18] for computing them in the cubic 
case. 


3. Uniform approximation. 


3.1 Fixed knots. Let a<xj)<x2<... <x, <b be pre- 
scribed and let $n, K(X, 0. X) be the class of splines of 


degree n with knots fxj}) defined in §2.1. 


Problem 3.1, Given fe C[a,b], find s*(t)e Sy , 
(x),.-+,%,) such that : 


Pd 
le- sll, pee I¢t) -s*yl <lle- sll, 


for every se Sane (X),--~)X,) ‘ 

The theory of this best approximation problem has been 
surveyed in [22] in this volume. In particular, we recall 
that a best approximation exists although in general it may not 
be unique. 

For computational purposes, we are forced to replace 
Problem 3.1 with the following discretized version. 


Problem 3.2. Let N bea positive integer and let 


D=fa+ih}N , where h =(b-a)/N. Given fe C[a, b] 
ind Ss” ¢ Sn, k(X),--.,%,) such that 
]£-s*ll, = max_ |f(x) - s*(x)| < lf - sll 
Dy” xe D ea D 


for every s€Sy x (X],.--, Xp) - 

The existence of best approximations in this sense is 
also easily verified, and most of the characterization theory 
quoted in [22] also carries over. Of special interest is the 
following de la Valleé Poussin type of result for splines (cf. 
Theorem 2. 8 of [22]). 





Theorem 3.3. Suppose there exist ies and 


{ante with xy>0 and a<tp<ty<... < theka <b 
Satisfying 
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(3.1) (f£-s)(t) =(-1)' &,, te: 2) Sembee (Bee 
Then 
min ys A, y(D) silp-sl, 
where 
A, (P) = min lle - sll, : 
SE Bn, K(X) +++) X) 


The importance of this result here, is that if some s 
exists satisfying (3.1) then upper and lower bounds for the 
deviation An, ,(D)_ are immediately available. For certain 
choices of te) pete the construction of such splines is 


quite easy. We calla set Mog = sa y} cD satisfying 
(0)- ,(0) (0) 
(3. 2) BS Sho RaSh eB 


admissible provided there exists a unique spline Sol t) € Sn, k 
(x),-.+,X%,) such that 


(0) dre aap) ~ : 
(3.3) aft) 4-0 Xo =H; ), Jal, 2,..., ntk+2, 


J 4 
for some \g #0; it is shown in [21] that a set M = (pt #2 


is admissible provided there exist {yp pte CM such that 
Yj <Xi <Yntie1, i =1,2,...,k. Now, as mentioned above, 
Sn, 7 64 ae »Xy) is a finite-dimensional linear space, and 


hence if {b, t) ii provide a basis for it, then we may 


m1 nt 
seek s satisfying (3.3) in the form s(t) = Lae aj b(t), 
i.e. , we solve the system of equations 
n+k+1 
0 at ih OD 
(3.4) 2 a b(t) #(-N dy =H’), 


$al,2,...,mene2, tor LaF eng ny. 

These ideas were exploited in [21] to derive a Remez 
type single-exchange method wherein M is adjusted so that 
larger and larger A's are obtained thus providing a sharp 
lower bound for An, K( D) SAn, k= min l f-s|l.. 

sé Sn, k( Xp + Xp) 
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An upper estimate for An, k is provided by lz = sllp+ o( h), 
where w is the modulus of continuity of f. 

Since at each step of this algorithm the linear system 
(3.4) must be solved, it is important to choose a good basis 
for Sp, k(x1,-+-+,X,) which avoids ill-conditioned matrices, 
It has been my experience that the basis exhibited in ( 2. 3) 
provides excellent results. Convergence properties of this 
algorithm have not been established, but it has been shown 
that the sequence of X's produced are increasing (see [21]). 

An alternate attack on Problem 3. 2 due to Esch and 
Eastman [9] consists in considering it as a linear program- 
ming problem. Standard computational techniques for such 
problems can then be applied. The linear programming method 
is recommended since the Remez algorithm discussed above 
is subject to restrictions which can cause it to terminate prior 
to obtaining the best approximation with respect to \|- | D 
(see [21]) . 

These algorithms have been tested on a number of 
smooth and unsmooth functions. For purposes of comparison 
we include an example in which we approximate the roof func- 


tion 
e do redse 
r(t))= 


l-t otherwise 
on [0,1] with N =100 and D = {i/100}10° ; 
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The following general remarks seem warranted by the 
experiments: (1) for smooth functions, splines in 8n,k do 
not do appreciably better than polynomials with the same num- 
ber of parameters (i.e., ntk+1); (2) for less smooth func - 
tions which are more or less piecewise smooth, splines are 
by far superior provided the knots are properly located. 


3.2 Free knots. In this section the class of approximating 
splines will be taken to be 


8n,k = {s( t) | there exist a =X <x, Siete <x) =) 
and integers M),+++,™M), with 1 sm, <nt+l and 

> m,= k such that s(t) € m,, in each of the intervals 
i=l n-m; 

(5,4) while sé C in an open neighborhood 
of x, Pat once ath ‘ 


This is the class of all splines of degree n with some k 
knots in the interval [a,b] allowing multiple knots up to 
multiplicity ntl (see [22], §§2, 3). 


Problem 3.4. Given fe Cla,b] find s*« Sy ~ such 
that Fi t 
f=." |= f(t) - < |lf - 
lf - s* ll 0 ere, 9 s"(t)| <lf- slo 


for every S€ Sp k- 

Theoretical results for this problem have been surveyed 
in [22]. As with Problem 3.1, in practice one must replace 
Problem 3. 4 with a discretized version as in Problem 3. 2. 
Research on effective methods for obtaining best spline ap- 
proximations with free knots is still in progress. Esch and 
Eastman [9, 10] have reported quite successful experience 
with a method based on solving the fixed knot problem coupled 
with a knot adjustment procedure. Their method is predicated 
on the assumption that there exists a best approximating 
S€8n,k such that f-s exhibits a normal alternating char- 


acter, i.e., there exist a<t) <... <th4p42 <b such that 
(£- s)(t) =(-1)°S By yy dal, 2,..., ntk42(h =), 
where B = min f=s : This, however 
ny k oe —_ Heo : 
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remains, a conjecture (cf. discussion in §3.2 of[22]). All 
of the examples treated in[9, 10] possessed alternating 
best approximations, but on the other hand, they alsosatisfied 
s(nt1)(t) > 0 which is known to be sufficient to guarantee 
the normal alternating characterization (see Theorem 5.2 and 
the discussion in §5.5 of [22]). 

Since the details of the Esch and Eastman algorithm 
are rather complicated, we do not include them here. We 
shall, however, quote a few of their examples and discuss 
some of their computational experience. 

The first example is approximation of the analytic func- 
tion f(x) =e* on [0,1] (see Table 1 below). It is seen that 
increasing the degree of the spline was more effective than 
adding knots. 

The second example involves approximating Nx on 
[0,1] (Table 2). Since Nx is not so well behaved near 0 , 
it is perhaps not surprising that the splines were more effec- 
tive than polynomials with a comparable number of parameters. 


4. Lz approximation. 


4.1 Fixed knots. For prescribed a <x] <X2 <<... Sx, <b 
let Sy x(x],-++,%,) be as in §§2.1, 3. 2. 


Problem 4.1, Given fe C[a,b] determine s*¢ Sn, k 
(x),...,x,) such that 


It -s*lla= JP (xt -s%(y]? at<lle- sll, 


for every se¢ Sn, KIX, »Xp) - 

Very little theory is available for best Lg approxima- 
tions (see §4 of [22], pp. 65-85, this volume). Since Sn k 
Chi sce »X,) has dimension n+k+l, the computational prob~ 
lem is a standard least squares problem, i.e., given a basis 
{b,(t) pte, determine an orthonormal basis {y,( ty} DtK+L 
e.g., by Gram-Schmidt. The choice of initial basis and ortho- 
normalization procedure are critical, however. A complete 
discussion as well as an explicit algorithm (and program) can 
be found in de Boor and Rice {7]. It has been tested exten- 


sively and has provided smooth spline approximations ina 
variety of problems. 
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k = Number of Knots 
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4,2 Free knots. Let Sn, k be the class of splines of degree 
n with some k knots (allowing multiplicities) defined in 
§3.2, 


Problem 4.2 Given fe Cla, b] determine s* ¢ $n, k 
such that 


le - stl, = [Pint - s(t]? ar< If- sl, 


for every S€ Sn k- 
: As was the case for the uniform norm (§3.2), this is 

a nonlinear approximation problem which, however, does admit 

of a solution; see [8] where de Boor and Rice have carried out 

experiments on an algorithm which finds a solution with fixed 

knots and then adjusts them and iterates. Details and a 


discussion of their experience can be found in [8]. 
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